2-2

Quantum Motion of a Particle in a Box

Nicholas Wheeler
Reed College Physics Department
13 August 2008

Introduction

It is by way of preparation for discussion of "quantum Liossajous figures" and "quantum vorticity" that I look here to aspects—specifically, dynamics aspects—of the simplest and most familiar of all quantum mechanical systems: the "particle-in-a-box" problem. I work first in one dimension, then by straightforward generalization in two dimensions.

One-dimensional theory

The simple essentials of the subject at hand are treated in §2.2 (pages 30–37) of Griffiths' Introduction to Quantum Mechanics (2nd edition, 2005). Some of what I will have to say parallels work reported in Cimarron Wortham's Reed College thesis: "Motion of launched wave packets in the infinite square well" (2004).

We are interested in properties of the solutions of

-((ℏ^2)/(2m))(∂_(x,x))ψ[x,t]=i ℏ (∂_t) ψ[x,t]

that conform to the boundary conditions

ψ[0,t]=ψ[a,t]    :    all  t

If we adopt units in which

ℏ=2m=a=1

we have

-(∂_(x,x))ψ[x,t]=i  (∂_t) ψ[x,t]
ψ[0,t]=ψ[1,t]    :    all  t

which is easier to work with. If we assume the wave function to have separated structure

ψ[x,t]=Ψ[x]*(e^(-i ω t))

we find ourselves looking for solutions of

-(∂_(x,x))Ψ[x]=ω Ψ[x]
 Ψ[0]= Ψ[1]=0

and are led immediately to the orthonormal eigenfunctions

Ψ[x_,n_]:=(Sqrt[2]) Sin[n π x]

Table[(∫_0^1)Ψ[x,m]Ψ[x,n]\[DifferentialD]x,{m,1,4},{n,1,4}]//MatrixForm

Plot[{Ψ[x,1],Ψ[x,2],Ψ[x,3],Ψ[x,4]}, {x,0,1},Ticks->{{1},{(Sqrt[2]),-(Sqrt[2])}},AxesLabel->{"x"}, LabelStyle->Directive[Bold], PlotLabel->"First four eigenfunctions"]
Plot[{((Ψ[x,1])^2),((Ψ[x,2])^2),((Ψ[x,3])^2),((Ψ[x,4])^2)}, {x,0,1},Ticks->{{1},{2,-2}},AxesLabel->{"x"}, LabelStyle->Directive[Bold], PlotLabel->"Associated Probability Densities"]

From

-(∂_(x,x))Ψ[x,n]==(n^2) (π^2)Ψ[x,n]

we are led to define the "eigenfrequencies"

ω[n_]:=(n^2)(π^2)

and to construct the "buzz factors"

BuzzFactor[t_,n_]:=(E^(-I ω[n] t))

The time-dependent eigenfunctions are give now by

ψ[x_,t_,n_]:=Ψ[x,n]*BuzzFactor[t,n]

which we verify:

-(∂_(x,x))ψ[x,t,n]==I (∂_t) ψ[x,t,n]

The following figures illustrate the motion of the complex-valued functions ψ[x, t, n], {n = 1, 2}:

Animate[ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,Ψ[x,1]Cos[-ω[1](t/(200π))],Ψ[x,1]Sin[-ω[1](t/(200π))]}}, {x,0,1}, PlotRange->{{0,1},{-2,2},{-2,2}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}},
PlotLabel->"Motion of (Subscript[ψ, 1])(x, t)"],{t,0,400}]

Animate[ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,Ψ[x,2]Cos[-ω[2](t/(200π))],Ψ[x,2]Sin[-ω[2](t/(200π))]}}, {x,0,1}, PlotRange->{{0,1},{-2,2},{-2,2}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}},
PlotLabel->"Motion of (Subscript[ψ, 2])(x, t)"],{t,0,400}]

The functions ψ[x, t, n] move (in the sense illustrated above), but the associated probability densities are static. Observable quantum motion arises when such functions are superimposed, as an interference effect. We look to a typical (normalized) superposition:

Motion of a simple superposition

ψ[x_,t_]:=(Sqrt[(2/3)])ψ[x,t,1]+(Sqrt[(1/3)])ψ[x,t,2]

(∫_0^1)((ψ[x,0])^2)\[DifferentialD]x

ComplexExpand[ψ[x,t]]

Realψ[x_,t_]:=((2 Cos[(π^2) t] Sin[π x])/(Sqrt[3]))+(Sqrt[(2/3)]) Cos[4 (π^2) t] Sin[2 π x]
Imagψ[x_,t_]:=-((2 Sin[(π^2) t] Sin[π x])/(Sqrt[3]))-(Sqrt[(2/3)]) Sin[4 (π^2) t] Sin[2 π x]

Animate[ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,Realψ[x,(t/(200π))],Imagψ[x,(t/(200π))]}}, {x,0,1}, PlotRange->{{0,1},{-2,2},{-2,2}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}},
PlotLabel->"Motion of a Simple Superposition"],{t,0,400}]

Animate[Plot[((Realψ[x,(t/(200π))])^2)+((Imagψ[x,(t/(200π))])^2),{ x,0,1},PlotLabel->"Sloshing of the Associated Probability Density"],{t,0,400}]

Motion of a more complex superposition

We acquire need at this point of the "beta distribution"

PDF[BetaDistribution[α,β],x]

which is defined for all real non-negative values of α and β, and clearly vanishes at both ends of the unit interval:

PDF[BetaDistribution[α,β],0]
PDF[BetaDistribution[α,β],1]

and in between manages for many (α, β) values to look quite "Gaussian":

Manipulate[Plot[PDF[BetaDistribution[α,β],x], {x,0,1}, PlotRange->{0,12},PlotLabel->"Parameter Dependence of the β Distribution"], {α,2,100},{β,2,100}]

The beta distribution provides a handy description of what might be called a "boxed Gaussian"—a normal distribution with its boundary-condition-violating tails cut off. (Wortham supplies a proof (pages 10–12 of his Reed thesis) that the beta distribution—despite its purely algebraic structure—becomes Gaussian in the small-variance limit.)

Let us look to the quantum motion that evolves from the initially real-valued wavefunction

Initialψ[x_]:=(Sqrt[PDF[BetaDistribution[50,100],x]])

Mean[BetaDistribution[50,100]]
(Sqrt[Variance[BetaDistribution[50,100]]])//N

Plot[{Initialψ[x],((Initialψ[x])^2)},{x,0,1}, PlotRange->{0,12}, Ticks->{{(1/3),1},{2,4,6,8,10,12}},PlotLabel->"Initial Wavepacket and Associated Probability Density"]

To launch such a wavepacket into motion, we write it as a weighted sum of eigenfunctions,  and then set each of those abuzz:

Coefficients=Table[{n,NIntegrate[Initialψ[x]*Ψ[x,n],{x,0,1}]},{n,1,40}]

ResolvedInitialψ[x_]:=(∑_(n=1)^31)Coefficients[[n]][[2]]Ψ[x,n]

The following figure

Plot[{Initialψ[x],ResolvedInitialψ[x]},{x,0,1}, PlotRange->{0,4},PlotLabel->"Compare at Time t = 0"]

provides convincing evidence that our  ResolvedInitialψ[x] does indeed provide an excellent approximation to Initialψ[x]. As also does the following proof of normalization:

NIntegrate[((ResolvedInitialψ[x])^2),{x,0,1}]

If we now turn on the buzz factors we are led to define

Evolvedψ[x_,t_]:=(∑_(n=1)^31)Coefficients[[n]][[2]]ψ[x,t,n]

The complex-valued function thus defined is fairly complicated:

ComplexExpand[Evolvedψ[x,t]]

We isolate the real and imaginary parts

RealEvolvedψ[x_,t_]:=ComplexExpand[Re[Evolvedψ[x,t]]]
ImagEvolvedψ[x_,t_]:=ComplexExpand[Im[Evolvedψ[x,t]]]

and use that information to construct a display of the motion of the complex-valued wavepacket:

Manipulate[ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,RealEvolvedψ[x,(t/(200π))],ImagEvolvedψ[x,(t/(200π))]}}, {x,0,1}, PlotRange->{{0,1},{-3,3},{-3,3}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}}],{t,0,400,1}]

The animation demonstrates, by the way, the recurrence that follows from the fact that

ω[n]==(n^2)ω[1]

The equation shows that the (n^th) mode completes (n^2) cycles in the time it takes the fundamental mode to complete one. Any wavepacket assembled from such buzzing eigenfunctions therefore will return periodically to its original form, the recurrence time being given by

(Subscript[T, recurrence])=((2π)/(ω[1]))

Here is a more "physical" derivation of the same result:

(Subscript[T, recurrence])=((2π)/(ω[1]))=((2π ℏ)/GroundstateEnergy)=2π ℏ((((((π^2)(ℏ^2))/(2m (a^2)))))^(-1))

Simplify[2π ℏ((((((π^2)(ℏ^2))/(2m (a^2)))))^(-1))]/.{ℏ->1, a->1,m->1/2}

All the complicated phase information—evident in the preceding figure—is lost when we construct the moving probability density (which itself is quite complicated, and again shows recurrence):

EvolvedProbabilityDensity[x_,t_]:=((RealEvolvedψ[x,t])^2)+((ImagEvolvedψ[x,t])^2)

Manipulate[Plot[Evaluate[EvolvedProbabilityDensity[x,(1/600)(Subscript[T, recurrence])n]],{x,0,1},PlotRange->{0,12},PlotLabel->"Motion of the Probability Density"],{n,0,600,1}]

Using that moving distribution to compute the moving expected position (this will take about a minute)

Evaluate[NIntegrate[Re[x*EvolvedProbabilityDensity[x,(1/600)(Subscript[T, recurrence])110]],{x,0,1}]]

ExpectedPositionData=Table[NIntegrate[Re[x*EvolvedProbabilityDensity[x,(1/600)(Subscript[T, recurrence])  n]],{x,0,1}],{n,0,600,1}];

ListPlot[ExpectedPositionData, Joined->True, PlotRange->{0,1},PlotLabel->"Motion of the Expected Position", AxesLabel->{"Time","Position"}]

Consider the enormous amount of information that has been discarded in the construction of the preceding figure—phase information which could never be reconstructed on the basis of the figure. 

We work in units in which the recurrence time is 600. It is curious, therefore, that the expected position executes (in this instance) three cycles during the span of a single recurrence cycle. The following superposition of figures casts some light on the situation

Plot[{EvolvedProbabilityDensity[x,(0/3) (2/π)],EvolvedProbabilityDensity[x,(1/3) (2/π)],EvolvedProbabilityDensity[x,(2/3) (2/π)]},{x,0,1},PlotRange->{0,12},PlotLabel->"Densities at Three Times, Superimposed"]

and the following figure casts an even brighter light:

Table[{"((Subscript[T, recursion])/3)"k,ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,RealEvolvedψ[x, (2/(3π))k],ImagEvolvedψ[x, (2/(3π))k]}}, {x,0,1}, PlotRange->{{0,1},{-3,3},{-3,3}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}}]},{k,0,2,1}]//TableForm

Launched wavepackets

We have been working thus far with a wavepacket that was initially at rest, in the sense that initially ⟨p⟩ = 0, as I attempt to demonstrate: recalling that we have set ℏ = 1, we have

InitiallyExpectedMomentum=-I NIntegrate[ResolvedInitialψ[x]D[ResolvedInitialψ[x],x],{x,0,1}];

InitiallyExpectedMomentum

…which is to say (if we disregard Mathematica's complaint): the expected value of the momentum is zero.

To "launch" such a wavepacket, one hits it with the complex multiplier

(e^((i/ℏ) p x))    :    p=m v

which in our units becomes

(E^((I/ℏ) p x))/.{ℏ->1, p->(1/2)v}

Suppose we were to set v = 300:

(E^((I v x)/2))/.v->300

Hitting a real-valued initial state with such a "launch factor" produces a complex-valued initial state, as illustrated below:

ParametricPlot3D[{{x,0,0},{0,x,0},{0,0,x},{x,ResolvedInitialψ[x]*Cos[150x],ResolvedInitialψ[x]*Sin[150x]}}, {x,0,1}, PlotRange->{{0,1},{-4,4},{-4,4}}, Boxed->False, Axes->False, PlotStyle->{{Thickness[.01]},{Blue},{Blue},{Red}},
PlotLabel->"Launched (∴ Complexified) Initial State"]

The initial expected value of the momentum has now the non-zero value

InitialLaunchedMomentum=-I  NIntegrate[ResolvedInitialψ[x]*(E^(-I 150 x))D[ResolvedInitialψ[x]*(E^(I 150 x)),x],{x,0,1}];

InitialLaunchedMomentum

InitialLaunchedψ[x_]:=(E^(I 150 x))*Initialψ[x]

Experiment (see below) shows that the coefficients that contribute most importantly to the resolution of a launched wavepacket have orders in  the range

(momentum/π)±20

In the present instance we have

(150/π)//N

so (to establish the point) we look at

NIntegrate[Ψ[x,28]*InitialLaunchedψ[x],{x,0,1}]
NIntegrate[Ψ[x,48]*InitialLaunchedψ[x],{x,0,1}]
NIntegrate[Ψ[x,68]*InitialLaunchedψ[x],{x,0,1}]

We construct a complete set of "launched coefficients"

LaunchedCoefficients=Table[NIntegrate[Ψ[x,n]*InitialLaunchedψ[x],{x,0,1}],{n,28,68}]

Length[LaunchedCoefficients]

which are now complex. We separate their real from their imaginary parts

RealPartsLaunchedCoefficients=Table[Re[LaunchedCoefficients[[k]]],{k,1,41}]
ImagPartsLaunchedCoefficients=Table[Im[LaunchedCoefficients[[k]]],{k,1,41}]

ListPlot[{RealPartsLaunchedCoefficients,ImagPartsLaunchedCoefficients}, PlotStyle->{PointSize[0.02],PointSize[0.02]}]

ListPlot[{RealPartsLaunchedCoefficients,ImagPartsLaunchedCoefficients}, Joined->True]

We use this information to assemble

MovingLaunchedψ[x_,t_]:=(∑_(n=1)^41)LaunchedCoefficients[[n]]ψ[x,t,27+n]

RealMovingLaunchedψ[x_,t_]:=ComplexExpand[Re[MovingLaunchedψ[x,t]]]
ImagMovingLaunchedψ[x_,t_]:=ComplexExpand[Im[MovingLaunchedψ[x,t]]]

LaunchedProbabilityDensity[x_,t_]:=((RealMovingLaunchedψ[x,t])^2)+((ImagMovingLaunchedψ[x,t])^2)

UnlaunchedInitial=Plot[Evaluate[EvolvedProbabilityDensity[x,0]],{x,0,1}, PlotRange->{0,12}];
LaunchedInitial=Plot[Evaluate[LaunchedProbabilityDensity[x,0]],{x,0,1}, PlotRange->{0,12}];

The initial launched distribution, though assembled now from a different set of (complexly weighted) eigenstates, again reproduces quite accurately the initial unlaunced distribution:

Show[UnlaunchedInitial,LaunchedInitial]

I show below the appearannce of the lunched probability density at each of 24 equally-spaced times during a single recurrence cycle:

Table[{k"((Subscript[T, recurrence])/24)",Plot[Evaluate[LaunchedProbabilityDensity[x,(2/π)(k/24)]],{x,0,1}, PlotRange->{0,12}]}, {k,0,24,1}]//TableForm

The following (fairly time-consuming) calculations evaluates the expected position at each of 6000 times during a single recurrence cycle:

ExpectedPositionDataAfterLaunch=Table[NIntegrate[Re[x*LaunchedProbabilityDensity[x,(1/600)(2/π)  n]],{x,0,1}],{n,0,600,1}];

ListPlot[ExpectedPositionDataAfterLaunch, PlotRange->{0,1}]

ListPlot[ExpectedPositionDataAfterLaunch, Joined->True,PlotRange->{0,1}]

One cannot help but notice the curious sub-periodicities, which I will not attempt to explain.

Quantum motion in the classical limit

Reinstating all physical parameters, we notice that the Schrödinger equation

-((ℏ^2)/(2m))(∂_(x,x))ψ[x,t]=i ℏ (∂_t) ψ[x,t]

can be written

-(∂_(x,x))ψ[x,t]=i ((2m)/ℏ) (∂_t) ψ[x,t]

which reduces to the equation with which we have been working if one adjusts the time-scale:

t→θ=(ℏ/(2m))t

One expects the quantum physics to look ever more classical as one approaches the limit  m ↑ ∞  (or equivalently: ℏ ↓ 0). Which amounts to shrinking the time scale

t→θ=α t   :   α≪ 1

We have encountered the familiar fact which lies at the base of the Feynman formalism—namely, that "quantum mechanics is briefly classical." I propose therefore to look to the short-time motion of a launched wavepacket (in sufficiently long times it will execute—exactly—all the complicated gyrations described above, and will, in particular, display recurrence).

Notice in this connection that in the most recent of the preceding figures the value of ⟨x⟩ oscillates back and forth during the brief interval

0<t<(20/600)(Subscript[T, recurrence])

with seemingly constant frequency and slowly diminishing amplitude. We take a magnified look at the quantum mechanics during that interval.

Table[{(k/20)"((Subscript[T, recurrence])/30)",Plot[Evaluate[LaunchedProbabilityDensity[x,(2/π)(k/600)]],{x,0,1}, PlotRange->{0,12}]}, {k,0,20,1}]//TableForm

ExpectedPositionDataShortTerm=Table[NIntegrate[Re[x*LaunchedProbabilityDensity[x,(1/6000)(2/π)  n]],{x,0,1}],{n,0,200,1}];

ListPlot[ExpectedPositionDataShortTerm, PlotRange->{0,1}]
ListPlot[ExpectedPositionDataShortTerm, PlotRange->{0,1}, Joined->True]

I draw attention once again to how little quantum information survives in the statement of this semi-classical result.

The figure shows the motion of the expected position during the first  (200/6000)=(1/30) of a recurrence cycle. The period of the oscillations appears to be about (1/63) of a cycle. That number can be understood as follows: in physical variables, one has

Quantum Recurrence Time=((2π)/(ℏ(Groundstate Energy)))=((4π m (a^2))/((π^(2 ))ℏ))

((Classical Period)/(Quantum Recurrence Time))=((2a)/v)(((((4π m (a^2))/((π^(2 ))ℏ))))^(-1))

Simplify[((2a)/v)(((((4π m (a^2))/((π^(2 ))ℏ))))^(-1))]

%/.{a->1,ℏ->1,m->(1/2)}

The computed data, on the other hand, supplies (by best informal estimate)

((Classical Period)/(Quantum Recurrence Time))=(63/6000)

Solve[(π/v)==(63/6000),v]//N

This agrees very well with the fact that at the outset we set v = 300.